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Abstract 

We present a detailed account of the technical aspects of stochastic quantum molecular dynamics, an ap- 
proach introduced recently by the authors [H. Appel and M. Di Ventra, Phys. Rev. B 80 212303 (2009)] to 
describe coupled electron-ion dynamics in open quantum systems. As example applications of the method we 
consider both finite systems with and without ionic motion, as well as describe its applicability to extended 
systems in the limit of classical ions. The latter formulation allows the study of important phenomena such 
as decoherence and energy relaxation in bulk systems and surfaces in the presence of time-dependent fields. 



1. Introduction 

Time-dependent density-functional theory 
(TDDFT) calculations are currently enjoying a 
large popularity due to their efficiency and success in 
describing low-lying excitation energies in molecular 
systems pQ. In addition, many applications have 
been investigated with TDDFT. Examples include 
electronic transport [3l HI [5j, nonlinear optical 
response [5J, or atoms and molecules in strong laser 
fields [7J[B]. In the latter cases, the time-dependent 
Kohn-Sham (TDKS) equations are usually evolved 
in real-time. However, the majority of these stud- 
ies pertains to the description of closed quantum 
systems, since the corresponding TDKS equations 
describe a set of N particles evolving coherently 
in time. Q On the other hand, most experimental 
situations involve some level of energy dissipation 
and/or decoherence induced by either some envi- 
ronments to which the given system is coupled, or 
the measurement apparatus itself which necessarily 
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projects non-unitarily the state of the system onto 
states of the observables. This is generally true for 
both electrons and ions, so that a first-principles 
description of their coupled dynamics in the presence 
of one or more environments is of fundamental 
importance in order to describe phenomena and 
compare with experiments. At this point, it is worth 
noting that present quantum molecular dynamics 
(QMD) approaches, (e.g., the Born-Oppenheimer, 
Ehrenfest or Car-Parrinello methods) either do not 
allow excited states dynamics (Born-Oppenheimer 
and Car-Parrinello methods) or, if they do (e.g., 
Ehrenfest QMD), they do not permit electronic 
coupling to external environments. Indeed, in all 
these approaches, energy dissipation and thermal 
coupling to the environment are usually described 
with additional thermostats coupled directly to the 
classical nuclear degrees of freedom, which fall short 
of describing the numerous physical phenomena 
associated with decoherence and energy dissipation. 

In order to overcome these shortcomings, we have 
recently introduced a novel time-dependent den- 
sity functional approach based on stochastic time- 
dependent Kohn-Sham equations [13], where we allow 
the coupling of both electrons and (in principle quan- 
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turn) ions with external baths. This approach - we 
have named stochastic quantum molecular dynamics 
(SQMD) - extends the previously introduced stochas- 
tic time-dependent-current density-functional theory 
(STDCDFT) HUD] to the coupled dynamics of elec- 
trons and ions. The latter was formulated to account 
for electrons interacting with external environments, 
without however including atomic motion. There- 
fore, SQMD combines and improves on the strengths 
of STDCDFT and present QMD methods by greatly 
expanding the physical range of applications of these 
methods. 

Clearly, from a practical point of view the present 
approach suffers - like all density-functional theory 
(DFT) based methods - from our limited knowledge 
of the properties of the exact exchange-correlation 
functional. Furthermore, in the present case, the ex- 
act functional depends not only on the electronic de- 
grees of freedom, but also on the ionic and bath(s) de- 
grees of freedom [T3]. Nevertheless, due to the weak 
system-bath(s) coupling assumption of the present 
theory, as well as the limited number of systems 
where quantum nuclear effects are of disproportion- 
ate importance, we may start by considering the limit 
of SQMD to classical nuclei and adopt the available 
functionals of standard closed-system TDDFT. Like 
in any other practical application of DFT, it is the 
predictions that we obtain and comparison with ex- 
periments that will be the ultimate judge of the range 
of validity of the approximate functionals used. 

In Ref. (T3] we have outlined the details of the proof 
of the theorem at the core of SQMD, and provided a 
simple example of the relaxation dynamics of a finite 
system (a molecule) prepared in some excited state 
and embedded in a thermal bath. However, there are 
some technical details behind an actual implementa- 
tion of this approach we have not reported yet, and 
which are nonetheless important if one is interested 
in using this method for practical computations. In 
this work we then present all the technical aspects 
for a practical implementation and use of SQMD. In 
addition, we present the theory behind its applicabil- 
ity to extended systems which is of great importance 
in the study of decoherence and energy relaxation in 
bulk systems and surfaces. We are in the process of 
implementing SQMD for extended systems and we 



will report these results in a forthcoming publication 

The paper is organized as follows. In Section [2] we 
give an introduction to the theory of stochastic quan- 
tum molecular dynamics. For completeness, this in- 
cludes general aspects of open quantum systems as 
well as the basic theorem of SQMD. In Section [3] we 
discuss the aspects of a practical implementation of 
SQMD. Finally, in Section [4] we illustrate with some 
examples the application of SQMD to finite systems 
with and without ionic motion, and outline its exten- 
sion to periodic systems. Conclusions are reported in 
Section [5] 

2. Theory 

2.1. Stochastic Schrddinger equation 

In the following, we consider an electron-ion many- 
body system coupled to a bosonic bath. For simplic- 
ity, we will consider only a single bath, but the for- 
mulation is trivially extended to the case of several 
environments. The total Hamiltonian of the entire 
system is then 



H = H s ®I B + Is®H B + \H SB - 



(1) 



Our system of interest is described by the many-body 
Hamiltonian Hs and the environment degrees of free- 
dom are given in terms of Hg- The interaction of the 
system with the environment is given by the Hamil- 
tonian Hsb and is assumed to be weak in the sense 
that a perturbation expansion in terms of this cou- 
pling can be performed. With A we denote the cor- 
responding coupling parameter for the system-bath 
interaction. 

The total system described by the Hamiltonian H 
follows a unitary time-evolution, which can be for- 
mulated for pure states either in terms of the time- 
dependent Schrodinger equation (TDSE), with h = 1 



= H(tMt), 



(2) 



or, alternatively for mixed states, in terms of the 
Liouville-von Neumann equation 



dt 



p{t) = -i H(t),p{t) 



(3) 
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where p is the statistical operator. 

Since H is a many-body Hamiltonian and we have 
to deal, in principle, with infinitely many degrees of 
freedom (due to the bath) it is not possible to solve 
Eq. ([2]) or ([3]) in practice, except for a few simple 
model cases. In addition, in most cases of interest the 
microscopic knowledge about the bath is limited and 
only its macroscopic thermodynamic properties are 
known, e.g., one typically assumes that the bath is in 
thermal equilibrium. However, we are only interested 
in the dynamics of the system degrees of freedom. It 
is therefore desirable to find an effective description 
for the system only. 

To accomplish this we may trace out the bath de- 
grees of freedom at the level of the statistical opera- 
tor, namely we perform the operation ps = Ttb{p}, 
where ps is called the reduced statistical operator of 
system S. It is worth pointing out here that this 
procedure does not generally lead to a closed equa- 
tion of motion for the reduced statistical operator and 
one needs further approximations. Depending on the 
approximations involved, one may arrive at an effec- 
tive quantum master equation for the reduced den- 
sity operator ps [13 EH [H] . As we will discuss later, 
this approach has some drawbacks when used within 
a density-functional formulation, both fundamental 
- in view of the theorems of DFT - and practical, 
since solving for the density matrix is computation- 
ally more demanding than solving directly for state 
vectors. 

We therefore take here a different route. Instead of 
working with a derived/composite quantity like the 
statistical operator, we summarize briefly how the 
bath degrees of freedom can be traced out directly 
at the level of the wavefunction. The derivation that 
follows has been reported elsewhere in the literature 
(see, e.g., Ref. [18 ). We repeat some steps here for 
completeness and to clarify our starting point. 

To this end let us consider the set of eigenfunctions 
{Xti(xb)} of the bath Hamiltonian 

H B Xn{xB) = £nXn(xB), (4) 

with xb the bath's coordinates (including possibly 
spin), and expand the total wavefunction of Eq. ^ 
in the complete set of orthonormal states formed by 



{Xn(x B )}, namely^] 

*(x S , X B ; t) = ^ ^(Xs; t)Xn(x B ), (5) 
n 

with <j> n (xs',t) some functions (not necessarily nor- 
malized) in the Hilbert space of the system S. 

In order to see that in the presence of a bath the 
functions (p n (xs',t) form a statistical ensemble de- 
scribing the properties of the subsystem S, let us 
proceed as follows. For a general observable O5 of 
the system S we find after simple algebra (and using 
the orthonormality of the bath states Xu(xb)) 

(y(xs,XB;t)\ds\y(xs,x B ;t)) = 

Y,(Mxs;t)\d s \<j>n(xs;t)). (6) 

n 

Let us now normalize the functions 4> n {xs] t) by writ- 
ing 

ip„(x s ,t) = <p n {x s ;t)/p n {t) (7) 

with 

Pn(t) = (<j>n(%S;t)\(f) n (xs;t)), (8) 

which, according to Eq. ([5]) , is nothing other than the 
probability for the bath to be in the state Xu{xb)- 

We can now define the following statistical opera- 
tor 

Ps = ^2pn(t)\^n(xs;t))(ip n (x s ;t)\ ^ 

and immediately recognize that the average ([6| can 
be re-written as 

(y(x s ,x B ;t)\d s \*(x s ,x B ;t)) = Tr{p s O s }, (10) 



2 At first sight this expansion might seem formally similar to 
the factorization used for the Born-Oppenheimer (BO) approx- 
imation. However, in the BO case the expansion coefficients 
depend on the dynamical variables of both subsystems, elec- 
trons and nuclei. This dependence originates from the fact that 
the electronic Hamiltonian in the Born-Oppenheimer approxi- 
mation has a parametric dependence on the nuclear degrees of 
freedom. In the present case we assume that the partitioning 
of the total Hamiltonian in Eq. |TJ is such that expansion 
becomes exact. 
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namely, due to the interaction with the bath, the 
system S is necessarily in a mixture of states de- 
fined by the macrostate {p n (t),tp n (xs',t)}. We thus 
expect that the equation of motion for the repre- 
sentative wave-function ip(t) of the subsystem S to 
be "naturally" stochastic, namely we expect to find 
an equation of motion that provides the macrostate 

{p n (t),1pn(xs;t)}. 

In order to show this we follow the Feshbach 
projection-operator method [HI (20) and define the 
following projection operators 

Pa :=I S ® |Xn>(Xn|, (H) 

Qn~is®Y, lx*Xxfcl , (12) 

where 1$ is the identity in the Hilbcrt space of the 
system. The rationale behind the choice of the above 
operators is to obtain the equation of motion of a 
representative coefficient <f> n (xs;t). 

By acting with these projection operators on the 
many-body TDSE for the combined system and bath 
in Eq. ^ we arrive at 

id t P n V(t) = P n HP n y(t) + P n HQ n V(t) (13) 



idtQn^(t) = Q n HQ n ^(t) + Q n HP n V(t) (14) 



right hand side is a memory term that is recording 
the history of the time evolution^] 

Note that, up to this point, we have made no 
approximations, i.e., the time evolution given by 



Eq. (151 is still fully coherent. However, the solution 



of Eq. (15 1 is very involved and, apart from model 



systems, not feasible in practice. Furthermore, a so- 
lution would require the initial conditions for all the 
microscopic degrees of freedom of the bath. These 
cannot all be determined simultaneously by a mea- 
surement. In practice, one rather has only knowl- 
edge about macroscopic thermodynamic properties 
of the bath, like temperature and pressure. It is 
therefore common to perform the following additional 
approximations which are motivated by the form of 
the system-bath interaction and the thermodynamic 
properties of the bath: (i) due to the assumed weak 
coupling between system and bath the source and 
memory terms are expanded up to second order in 
the system-bath coupling parameter A, (ii) the bath 
and subsystem S are assumed to be uncorrelated at 
the initial time, (iii) a random phase approximation 
is performed for the phases in the source and memory 
tcrm^ and (iv) it is assumed that the bath degrees 
of freedom form a dense energy spectrum and are in 
local thermal equilibrium characterized by 



Pb = 



Tr(e 



1 

-0He 



(16) 



Equation (14) can be formally solved. Inserting the where /3 = l/fc^T. 



result back into Eq. ( 13 ) we obtain 



id t P9(t) = PHPP^(t) + PHQe^ QHQt Q^{0) 
i. I PHQe^^-^QHPP^iT) dr, 



o 



(15) 



where we have omitted the index n for brevity. The 



first term on the right hand side of Eq. ( 15 ) con 



tains only projections on the system manifold, and 
describes the coherent evolution of the system degrees 
of freedom. The second term is a source term that 
carries a dependence on initial conditions (Q^(0) are 
the initial conditions of all system's states except the 
one we are considering), and the third term on the 



3 These equations have a formal similarity to the quantum 
transport formulation introduced by Kurth et. al. 12 1 1 . How- 
ever, in this case the projection operators project on the real 
space regions of leads and central molecular device. Also the 
bath is fermionic in the quantum transport case (leads) and 
electrons can be exchanged between system and "bath" . This 
is in contrast to the present case where we consider bosonic 
baths and only energy and momentum can be exchanged be- 
tween system and bath. 

4 The random phase approximation invoked in the deriva- 
tion of the Markovian stochastic Schrodinger equation might 
seem at first sight surprising. The derivation of the Lindblad 
equation from the Markovian stochastic Schrodinger equation 
on the other hand shows, that both describe the exact same 
dynamics if the Hamiltonian does not depend on internal de- 
grees of freedom or any time-dependent or stochastic field (see 
Sec.|2~2l. 
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Let us then write the interaction Hamiltonian as 
HsB = '$2S a ®B a) (17) 

a 

where S a and B a are - in the most general case - 
many-body operators that act on the Hilbert spaces 
of the system and bath, respectively. In the following 
we will also assume that the average of the opera- 
tors B a vanishes on the n-th eigenstate of the bath, 
namely 

~]S a (Xn\B a \ X n) =0. (18) 



If this is not the case we simply redefine the system 
Hamiltonian via 



Xn\B a \ Xr, 



(19) 



and the interaction Hamiltonian as H' SB = Hsb — 

AE a S a (Xn \B a \Xn)- The term (xn\B a \ x n ) thus 
contributes to the unitary evolution of the system by 
renormalizing its eigenvalues (a typical example of 
this is the Lamb shift [TIjI I17j ). 

With these approximations in place, the source 
term can be regarded as a stochastic driving term. 
This is because, the system's state we have singled 



out in Eq. (11) now interacts with a (practically in- 



finite) large set of bath states densely distributed in 
energy. The previously coherent equation Eq. (15) 



then has to be regarded as a non-Markovian stochas- 
tic Schrodinger equation for the general state vector 

1>(t) = Mxsrt/iMxslWnixslt)) [IS] 

idMt) = H s1 p(t) + \ J2Ut)s a ip(t) 
t a 

-i\ 2 Y,l C a p{t-T)Sle- ift ^-^S^{T)dT 



a/3 



0(A 3 ), 



(20) 



where l a (t) are stochastic processes with zero ensem- 
ble average, l a (t) — 0, and correlation functions 



la(t)W) = 0, 



l* a {t)l {t') = C a0 (t-t'). 



(21) 



(22) 



Equation ( 20 1 is a general non-Markovian stochastic 



Schrodinger equation. Indeed, it still contains a time- 
integral over the past dynamics which is originating 
from the memory term of Eq. (15). Even though 



the theorem of SQMD could be formulated with non- 
Markovian baths we will focus in the following only 
on the Markovian limit 



C aP (t-l/)<xS a p5(t-l/), 



(23) 



namely, we consider baths that are ^-correlated. 
Physically, this means that the bath does not retain 
memory of the interaction with the system which is 
valid when the typical thermalization time-scales in- 
side the bath are much faster than the thermalization 
time-scales of the system. This approximation is well 
justified for a large number of bath degrees of free- 
dom. If this assumption does not hold, one has to 



resort to the solution of the more involved Eq. ( 20 ) 



By inserting the Markov approximation, Eq. (23), 
into Eq. (20) we then arrive at the stochastic 



Schrodinger equation in the Born-Markov limit 



idtfi>(t) =H S (t)m - l ~Y. sis a m 

a 

+ X)la(t)&V(*). 



(24) 



where the parameter A has been absorbed in the 
operators S a . The first term on the right hand 
side of Eq. (24) is the usual unitary evolution of 



the system under the action of the system Hamil- 
tonian Hg, the second term describes the dissipa- 
tion effects introduced by the bath and would in- 
deed make the probability density generated by ifi(t) 
decay in time. The last term, however, introduces 
fluctuations so that the norm of the state vector 
ip(t) averaged over the ensemble is conserved, namely 
(#P)) = l + 0(4 

Due to the stochastic nature of this equation, the 
stochastic process described by Eq. (24) has to be 



simulated in terms of an ensemble of state vectors 
ip(t). Each member ip(t) of the ensemble evolves dif- 
ferently in time due to the random variables l a (t) in 
the third term on the rhs. of Eq. ( 24 ). If we consider 
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an initial mixed state 

/3s(0)=5>£|^(0)>(V>*(0)|, (25) 
k 

where pi are the probabilities (with J^kPk = -0 °f 
finding the state V'fc(O) m the ensemble, the statistical 
average over all members of the ensemble allows us to 
construct the reduced density operator for the system 
degrees of freedom 



(26) 



Here, we use the symbol 7TT to indicate the statistical 
average over all members of the ensemble of state 
vectors ipk(t), namely the ensemble {ipKt)} of state 
vectors with initial conditions "0fc(O)- 

The expectation value of a general physical observ- 
able of the system S, 0$, can then be computed as 



in Eq. (10 1, i.e 



(Os) 



= Tr{6 s p s (t)) 

= Tr(dsY,P°k\Mt))(Mt)\) 



(27) 



= J2pl(Mt)\o s \Mt)), 



where the last step shows that the construction of 
ps(t) is not actually required: we can compute expec- 
tation values of observables directly from the wave- 
functions in the usual way, followed by a statistical 
average over all members of the ensemble of state vec- 
tors. It is also important to note that this approach 
provides directly the full distribution of the given ob- 
servable at any given time, provided we can compute 
a large enough set of realizations of the stochastic 
processes l a (t) (for an example of this see, e.g., Ref. 
|13j). From this distribution we can then compute all 
higher moments and/or cumulants (e.g., the variance, 
skewness, etc.) some of which are directly accessible 
exp eriment ally. 

2.2. Derivation of the Lindblad equation and stochas- 
tic Hamiltonians 

For many-body Hamiltonians which are not 
stochastic, namely they do not depend on internal 



degrees of freedom - Hamiltonians Hs ^ -ffsKIV'fe)}] 
- or do not depend explicitly on some stochastic field, 
like e.g., a stochastic thermostat [22] . it is possible to 
derive the Lindblad equation [H [Ml OH HI] from the 
stochastic Schrodinger equation (24 1. 



For notational clarity, let us denote in the following 
discussion with \ip) a single member of the stochastic 
ensemble {iV'fc)}- If we consider for simplicity the case 
of a single bath operator in Eq. ( 24 ) , and observe that 
in the Markovian limit 



W{t) 



l(t')dt' 



(28) 



is a Wiener process [15] with properties dW — 
and dW^dW = dt, we can formulate the stochastic 



Schrodinger equation ( 24 ) for a single bath in differ- 



ential form according to 



d\i>) = 



dt-iS\ip)dW. (29) 



Next, we employ Ito stochastic calculus in order to 
compute the following differential 

Unlike in normal calculus, we also have to keep the 
third term in the product rule above. This becomes 
necessary, since a statistical average over the Wiener 
increment dW^dW is proportional to dt, which will 
cause terms quadratic in dW to contribute to first 



order in dt. Inserting Eq. (29) and its Hermitian 



conjugate into Eq. (30) we arrive after elementary 
algebra at 



d\ip)(ip\ = - iS\ip)(ip\dW + h.c. 



dt- 



dt 



-i\H s ,\i>)(^\ 

+ s\^)Oj\sUwUw 

+ S\il>){ip\HsdWdt + h.c. 
+ l -S\i)){4>\S Ji SdWdt + h.c. 



H s \^){ijj\H s dt 2 + ^S\^)(^\S^Sdt 2 



(31) 



{h s , 



dt 2 . 
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In order to construct the statistical operator from 
the state vectors of the statistical ensemble {|V>j[.}}, 
we perform in the next step the statistical average 
over all members in the ensemble, i.e. 



dp = d\iP)(i>\. 



(32) 



Taking the properties dW = 0, dWdt = and 
dW^dW = dt of the stochastic process l(t) into ac- 
count, we see that only the second and third line in 



Eq. (31 ) contribute to first order in dt and we arrive 
at 



dp = 



6 8 ,\4)M]dt-±{&S,\il>)(1>\}dt 



We have thus shown that the stochastic 
Schrodinger equation of Eq. (24) and the mas- 
ter equation (34 1 lead to the same statistical 
operator, if and only if the Hamiltonian is not 
stochastic. However, in order to prove any DFT 
theorem one is led to consider the dynamics of the 
actual many-body system and that of any auxiliary 
one (including the Kohn-Sham system) with different 
interaction potentials, but reproducing the exact 
many-body density or current density. It is then 
at this stage that a choice has to be made - in 
the case of a many-body system open to one or 
more environments - regarding the basic equation of 
motion to work with. If we choose to work with a 



+ S\i;)(^\SUt + 0(dt 2 ). 



(33) 



quantum master equation of the type (34), then we 



At this point, note that this equation of motion is 
not necessarily closed for p = \tp)(ip\ because the 
first term on the right hand side of Eq. ( 33 ) is not 

equal to the commutator —i 



H, 



s,Ps 



unless Hs 

Hs[{\ip 3 k )}], or Hs does not depend on any stochas- 
tic field, or the system is in a pure state at all times 
- which would amount to the case S = 00 How- 
ever, if the Hamiltonian is stochastic, one has to deal 
with an ensemble of Hamiltonians, and the statisti- 
cal average of the first term on the right hand side of 



are assuming from the outset that the Kohn-Sham 
Hamiltonian is not stochastic. But this is an hypoth- 
esis that constitutes part of the final thesis, namely 
we have to prove that this statement is correct, not 
assume it a priori |27j . This issue does not arise with 
the stochastic Schrodinger equation ( 24 ) , because in 



that case we can consider all possible Hamiltonians, 
including those that are stochastic. 

In addition to the above important point, we also 
recall that for arbitr ary time-dependent operators 
S(t) and H s (t), Eq. (34) may not yield a positive- 



Eq. ( 33 ) involves also a statistical average over these 



Hamiltonians (see, e.g., Refs. [231 121)] ). 

For the moment being, let us assume that Hs 
•ffs'[{|V , fe)}] an< 4 furthermore that the Hamiltonian 
Hs does not depend on some external stochastic field. 
In this case we find 



dtps 



Hs,ps 



which is the well-known quantum master equation in 
Born-Markov limit (or Lindblad equation if the bath 
operators and the Hamiltonian, do not depend on 
time) [H [Ml US El- 



definite statistical operator at all times (see, e.g., 
Ref. [25]). This is a major limitation in practical 
calculations, since loss of positivity (which precludes 
a statistical interpretation of physical observables) 
should then be checked at every instant of time. Note 
that such a limitation docs not pertain to the stochas- 
tic Schrodinger equation which can be equally ap- 
f . plied to arbitrary time-dependent operators without 

\&S,ps\+Sps& (34) 

possible loss of positivity of the ensuing statistical 
operator. Therefore, the above two issues make the 
equation of motion of the statistical operator a less 
solid starting point for a DFT theory of open quan- 
tum systems. 



5 A further complication would arise if the operators 5 de- 
pended on internal degrees of freedom, i.e., S = ■S[{|V'fc}}]- in 
that case, the average over the statistical ensemble in the sec- 
ond and third terms on the right hand side of Eq. ( | 33 [ > has to 
be performed over the operators S as well. 



2.3. Theorem of Stochastic QMD 

We are now in a position to state the basic theorem 
of SQMD. Before doing this let us define the basic 
quantities we work with. The many-body system we 
are interested in consists of N e electrons with coor- 
dinates r = {rj} and N n — J2 s N s ,n nuclei, where 
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each nuclear species s comprises N s<n particles with 
charges Z s j, masses M s j, j — l...N SiTl , and coor- 
dinates R = {R SJ }, respectively. Their dynamics - 
subject to an arbitrary classical electromagnetic field, 
whose vector potential is A(i) - is described by the 
Hamiltonian 

H s {t) = f e (r,t) + W cc (r) + Uext.efei) + W en {r,R) 
+ T n (R, t) + W nn (R) + L> cxt ,„(R, t), (35) 

where T e (t) and T n (t) are the kinetic energies of 
electrons and ions, with velocities Vk(t) = [f>k + 
eA(? fc) t)]/m and V a {t) = [P a - Z a A(R a ,t)]/M a , 
respectively and f7 e xt,e(r> t), U ex t,n(R,t) the external 
potentials acting on electrons and ions. The particle- 
particle interactions are given by 



W ee (r) 
^nn(R) 



1 N Q 



47re 'rf lr 

— y 



j - r fei 



j<k 



h), 



Rfi 



W en (r,R) 



N n 

w nn (R„ - R/3 

a</3 

y^ y^ w cn (f fc - R a ). 

fe=l a=l 



We then define the charge current operator 

e 

2m 



y> fc (t),J(r-f k )} 



(36) 
(37) 



for the particle number, and 
J(x,i) =j(r,*) + 



J Q (R, <) 



(40) 



for the current density. To simplify the notation we 
have also denoted with x = {R, r} the combined set 
of electronic and nuclear coordinates, and we use the 
combined index a = {s,j} for the nuclear species. 

We now formulate the theorem for a single bath 
operator. It trivially extends to many operators. 
Theorem. — For a given bath operator S, many-body 
initial state ^(x, t = 0) (not necessarily pure) and ex- 
ternal vector potential A(x, t), the dynamics of the 
stochastic Schrodinger equation in Eq. ( 24 ) generates 



ensemble- averaged total particle and current densi- 
ties N(x,t) and J(x, t). Under reasonable physical 
assumptions, any other vector potential A'(x, t) (but 
same initial state and bath operator) that leads to 
the same ensemble-averaged total particle and cur- 
rent density, has to coincide, up to a gauge transfor- 
mation, with A(x, t). 

A sketch of the proof of this theorem can be found 
in the original paper (T3]. We thus refer the reader 
to this publication for more details. Here, we just 
mention an important point. As in Ref. [UJ HOj we 
are implicitly assuming that given an initial condi- 
tion, bath operator, and ensemble-averaged current 
density, a unique ensemble-averaged density can be 
obtained from its equation of motion: 



<9N(x, t ) 



dt 



V • J(x,f) 



(41) 



S^nS — -S^Sn — -nS^S 
2 2 



If we write this equation in the compact form 

cW(x, i) = -V • J(x, i) + Jb(x, t) (42) 

the above amounts to saying that Tb (x, t) is a func- 
z , „ , tional of N(x, t) and J(x, t), or better of J(x, t) alone, 

J a (R, t) = —— 2^ {Vp(t),d(R - Rp)}, (38) and that Eq p a d m its a unique physical solution. 



for the electrons, and 
Z, 



2M a 

Z a =Z l3 , P M a =M ll 

for the ions. The total particle and current density 
operators of the system can then be written as 



(39) 



Therefore, unlike what has been recently argued [12] . 
the density is not independent of the current den- 
sity, and our theorem establishes a one-to-one corre- 
spondence between current density and vector poten- 
tial. If this were not the case, namely that the parti- 
cle and current densities were independent functions, 
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then Js(x, t) would not be completely determined 
by the sole knowledge of N(x,t) and J(x, t) [57] . 

2.4. The limit of classical nuclei and Kohn-Sham 
scheme 

At this stage we may formulate a Kohn-Sham 
(KS) scheme of SQMD where an exchange-correlation 
(xc) vector potential A xc - functional of the initial 
states, bath operator(s), and ensemble-averaged cur- 
rent density - acting on non-interacting species, al- 
lows to reproduce the exact ensemble-averaged den- 
sity and current densities of the original interacting 
many-body system. The ensuing charge and current 
densities would contain all possible correlations in the 
system - if the exact functional were known. 

However, in the present case, we could construct 
several schemes - based on corresponding theorems - 
by defining different densities and current densities. 
For instance, we could collect all nuclear densities 
into one quantity as done in Ref. [29] . This, by no 
means is a limitation of this approach. Rather, it al- 
lows us to "specialize" the given schemes to specific 
physical problems. Instead, a much more serious lim- 
itation relates to the construction of xc functionals 
for the chosen scheme. Therefore, as anticipated in 
the introduction, we will restrict ourselves here to the 
limit of classical nuclei. 

Let us then assume that we know the vector po- 
tential A c ff that generates the exact current density 
in the non-interacting system. By construction, the 
system follows the dynamics induced by the stochas- 
tic Schrodinger equation (for a single bath operator) 



d\*Ks) = \-iH KS - l&s) \9 KS )dt , 

V 2 / (43) 

- iS\^ KS )dW 

where \^ks) is a Slater determinant of single-particle 
wave-functions and 

iv 2 
fj \- [gfc +eA cfF (f fc ,R,f)] 



k=l 



2m 



where A ext is the vector potential applied to the true 
many-body system, and Ah xc is the vector potential 
whose only scope is to mimic the correct dynamics of 
the ensemble-averaged current density, and we have 
lumped in it also the Hartree interaction potential in 
addition to the xc one. All these potentials depend 
on the instantaneous classical nuclear coordinates R. 

We immediately note that for a general many-body 
bath operator acting on many-body states one can- 



not reduce Eq. (43) to a set of independent single- 



is the Hamiltonian of non-interacting particles with 
A e ff(rfc,R,i) = A ext (f fc ,R,t) + A hX c(rfc,R,t), (45) 



particle equations, as done in the usual DFT schemes 
for closed systems. In other words, this would gen- 
erally require the solution of an equation of motion 
of Slater determinants, which is still computation- 
ally quite demanding. To see this point, suppose we 
have N particles and retain M single-particle states. 
We then need to solve for Cjjf — 1 elements of the 
state vector (with C% = M\/N\(M - N)\ and the 
— 1 comes from the normalization condition). In ad- 
dition, one has to average over an amount, call it m, 
of different realizations of the stochastic process. The 
problem thus scales exponentially with the number of 
particles. If this seems prohibitive let us also recall 
that a density-matrix formalism would be even more 
computationally demanding, requiring the solution of 
{C^J + 2) x (Cff - 1)/2 coupled differential equations, 
even after taking into account the constraints of her- 
miticity and unit trace of the density matrix. 

It was recently suggested in Ref. [30] that for op- 
erators of the type O = J2j Oj > namely operators 
that can be written as sum over single-particle opera- 
tors (like the density or current density), the expecta- 
tion value of O over a many-particle non-interacting 
state with dissipation can be approximated as a sum 
of single-particle expectation values of Oj over an 
ensemble of N single-particle systems with specific 
single-particle dissipation operators. In particular, 
it was found that the approximate single-particle 
scheme provides an excellent approximation for the 
current density compared to the exact many-body 
calculation. [30] We refer the reader to Ref. [30] for 
the numerical demonstration of this scheme and its 
analytical justification. 

From now on, for numerical convenience, we will 
then adopt the same ansatz which in the present case 
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reads, 



(^KS\0\^KS) 



KSI 



(46) 



with I0 3 



KS, 



? K S ) 



single-particle KS states solutions of 



.(p + eA eff (r,R,i)) z 



2 m 



dt 



(47) 



iSi p \^ KS )dW(t), 



with Sip an operator acting on single particle states. 



2.5. Model for the system-bath interaction 



The single-particle operators S 3 sp in Eq. (471 that 
we employ in the present work are given by the fol- 
lowing time-independent projectors^] 



x |Vr(E))(^ b (r)|, 



(48) 



where f D (e k ) 



1 + exp 



denotes the 



usual Fermi-Dirac distribution and <5fej (1 — 5kk') de- 
note Pauli blocking factors. The projectors in Eq. 



(48) cause a relaxation of the system back to the 
ground-state orbitals | (r) ) with a rate given by 
the rate constants 7jfc(r) (generally space and orbital 
dependent), while the temperature in the Fermi fac- 
tor is modeling the temperature of the bath. 

For the present paper, we further assume that the 
rates are independent of space and orbital indexes, 
i.e, 7jfc(r) oc 1/r, where r is a relaxation time. The 
bath operators of this model are sufficient for the 
illustration purposes of the present work but they 
clearly provide only a simplified picture of the full 
system-bath interaction. We emphasize here, that 
a rigorous form of the bath operators and the asso- 
ciated relaxation rates can always be derived from 
the microscopic form of the complete Hamiltonian 



of system and bath and their mutual interaction, 
Eq. ([I]). For example, in the case of a phonon bath, 
the system-bath interaction Hamiltonian could be 
taken into account in terms of e.g., a Frohlich in- 
teraction. In that case the relaxation rates can be 
extracted from the electron-phonon coupling matrix 
elements. The situation is much simpler is the case of 
a photon bath, where the Einstein rates of stimulated 
and spontaneous emmission can be used. 

2.6. Forces on ions 

Once we have the single-particle KS states and the 
corresponding Slater determinant $jfs(x, t) at hand 
we can compute the forces on the nuclei as |13j 



F a (*) 



(*jcs(x, t)\V a H KS (x, t)|*xs(x, t)} 



(49) 

for each realization of the stochastic process, [j Note, 
that this force is stochastic in nature since the wave- 
functions in the above expectation value are solu- 
tions of a stochastic Schrodinger equation. Since ap- 
proximations to the xc functional of the Kohn-Sham 
Hamiltonian may make the latter stochastic then the 
force one would obtain using a density matrix ap- 
proach - e.g., by solving the quantum master equa- 
tion (34) - would not be necessarily equal to the av- 



erage force obtained from Eq. ( 49 ) by averaging over 
the ensemble of realizations. 



3. Simulation Algorithms 

We have now outlined the general theory behind 
SQMD and we are ready to move on to the descrip- 
tion of its actual implementation. 

3.1. Real-time propagation 

The standard real-time propagation of the Kohn- 
Sham orbitals for a closed quantum system is based 



3 Cf. also to the examples in Refs. 30 31, 10, 13 . 



7 Note that Eq. j49| l is not the expression for the force 
one would obtain from the Hellmann-Fcynman theorem. This 
is because we are considering a system out of equilibrium. 
Rather, Eq. |49| is the total time derivative of the average 
of the ion momentum operator over the state of the system 
(see, e.g., Ref. [32]). 
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on numerical approximations for the time-ordered 
evolution operator 

/ rt+At \ 

U{t + At,t) =fexp ( -i J H K s(T)dT\. (50) 

There are several approaches employed in standard 
TDDFT computer packages, like, e.g., octopus [33J 



to evaluate Eq. (50) numerically. Here, we have cho- 



sen the Magnus propagator as basic building block 
for our stochastic simulations. [34] The Magnus se- 
ries |34j provides an exact expression for the time- 



evolution operator ( 50 1 as a time-unordered exponen- 



tial of so called Magnus operators in the form 



U(t + At, t) = exp ( Oi + Cl 2 + Cl 3 



(51) 



where the ttj are given in terms of time-integrals over 
nested commutators of the Hamiltonian at different 
points in time 



/t+At 
H K s(r)dT 

/•t+At /-ri 

&2=J J [HKs(Ti),H KS (T2)}dr 2 dn 



(52) 



3.2. Quantum Jump Algorithm 

We are now left with the actual solution of the 
stochastic Schrodinger equation (24). In past work, 



this has been done by directly integrating this equa- 
tion with standard approaches - e.g., with appropri- 
ately modified Runge-Kutta methods (see e.g. [T01I35] 
and references therein). These approaches are rea- 
sonable when we deal with a small number of ac- 
cessible states or short propagation times. However, 
they become increasingly unstable with an increasing 
number of states or for very long timescales, which 
is the case for realistic systems, like molecular struc- 
tures, surfaces or solids. As an alternative, we have 
thus adopted the quantum jump algorithm pioneered 
in the work of Diosi [35], Dalibard [37], Zoller and 
Gardiner and collaborators (38] [39] [40] as well as 
Carmichael |41J. At the price of introducing the prop- 
agation of auxiliary states, the quantum jump algo- 
rithm provides improved stability for systems with 
a large number of states/particles and, due to the 
piecewise deterministic evolution, also a stable prop- 
agation scheme for long timescales. 

This algorithm works as follows. Consider the 
deterministic time-evolution given by the follwing 
norm-preserving non-linear Schrodinger equation 



The time-integrals can be evaluated numerically with 
e.g., a Gauss-Legendre quadrature. In the simplest 
case, which is accurate up to second order in the time- 
step, one arrives at the exponential midpoint rule 



(53) 



U(t + At, t) = exp [9,1 j + O(A^) 

fix = -iH K s{t + At/2) + 0(At 3 ). 



Higher-order approximations can be easily derived 
from the Magnus series and appropriate quadrature 
points and weights, but experience shows that the 
second order gives a good balance between speed and 
accuracy for many applications. In the present work 
we use this approximation for the piecewise deter- 
ministic evolution that we are going to introduce in 
the next section. 



where the non-Hermitian Hamiltonian Us is given by 



Us = H s - -&S. 



(55) 



As before, the operators Hs and S denote the Her- 
mitian system Hamiltonian and the bath operator 
of Eq. §24$ The main objective of the quantum 
jump algorithm is to sample the stochastic process 
given by Eq. ( |24| in terms of a piecewise determinis- 
tic evolution, i.e. a set of deterministic time intervals 



generated by the evolution of Eq. ( 54 ) and action of 



the bath operator S between two consecutive time 
intervals. A central ingredient of the algorithm is 



8 For convenience we consider here the case of a single bath. 
A generalization to many baths is straightforward. 
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Figure 1: The figure illustrates the time-evolution as generated by the quantum jump algorithm. The lower track represents 
the piecewise deterministic propagation of the physical state which is intercepted at instances in time where the bath operator 
S acts on the state. The points in time where this takes place are determined by sampling a waiting-time distribution. The 
sampling is performed by propagating an auxiliary state (represented in the upper track) with a non-Hcrmitian Hamiltonian. 
Uniformly distributed random numbers are drawn and once the norm of the auxiliary state drops below the current random 
number the propagation of the physical and the auxiliary state is suspended. At this point in time the action of the bath 
operator on the physical state results in a new state which is then also used to initialize the auxiliary state for the evolution. 
The simulation of both states is then resumed again. 



a waiting-time distribution which determines when 
the jumps (i.e., actions of the bath operator) appear 
throughout the simulation. 

In order to sample the unknown waiting-time dis- 
tribution, an auxiliary set of wavefunctions (/>^ ux is 
introduced. The wavefunctions <^ ux are propagated 
with the non-Hermitian Hamiltonian alongside 
the actual states ipj. Since the auxiliary system 
evolves with a non-Hermitian Hamiltonian, the norm 
of the states (/>| ux is not preserved. It can be shown 
|16j that the decay of the norm of the auxilary wave- 
functions is related to the waiting-time distribution. 
The algorithm makes use of this fact and directly 
samples the waiting-time distribution on the fly from 
the norm decay. 

In terms of the Kohn-Sham system, the steps of 
the algorithm can then be summarized as follows 

1) Draw a uniform random number r/k <= [0, 1] for 
the Kohn-Sham Slater determinant 



2) Propagate TV auxiliary orbitals <j>j l 
non-Hcrmitian dynamics 



aux 



KS 



/aux 
3 ' 



under the 



,N 



3) Propagate the orbitals ipf S i 3 = 1 ■ ■ • N of the 
Kohn-Sham system with a norm-conserving dy- 
namics according to 



KS 



H K s-^^S + z\\S^f s \ 



4> 



KS 



4) If the norm of the Kohn-Sham Slater determi- 
nant drops below the drawn random number rjj , 
act with the bath operator(s) on the Kohn-Sham 
orbitals and update the auxiliary orbitals 



||Det{^ ux (i fe )}||< % ^ 



?Pf b (tk) 



5) Go to step 1) 

The piecewise deterministic evolution that is gen- 
erated by the steps of this algorithm is illustrated 
schematically in Fig. [T] 

Averaging at any given time over an ensemble of 
stochastic realizations allows then to obtain mean 
values of any physical observable. It is also impor- 
tant to realize that we have a full statistical ensemble 
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at hand. This allows to compute distributions of ob- 
servables, higher order moments, cumulants, etc. In 
this sense the ensemble of stochastic realizations gen- 
erated by the stochastic Schrodinger equation carries 
more information than the statistical operator, since 
the latter is merely a first order moment and higher 
order moments and cumulants cannot be computed 
easily from the first order moment (if at all) . 

We also emphasize here, that the interpretation of 
a single stochastic trajectory is not meaningful: the 
stochastic realizations have to be considered always 
as an ensemble. When averages over the stochastic 
ensemble are performed, the "convergence" of all ob- 
servables of interest has to be checked carefully by 
increasing the number of realizations of the stochas- 
tic process. 

Note that without further constraints the action of 
the bath operator in step 4) of the algorithm can in 
principle lead to a loss of orthogonality. For exam- 
ple, all orbitals of the Slater determinant could relax 
to the same orbital shape. The system could loose in 
this way its fermionic character. In order to maintain 
the fermionic nature of the Kohn-Sham state vector, 
we have to ensure that the orbitals of the Kohn-Sham 
Slater determinant remain orthogonal. To achieve 
this, we perform an orthogonalization of the orbitals 
after each action of the bath operators. This orthog- 
onalization can be thought of as being part of the 
definition of the action of the operator S. 

From our numerical experience so far, the waiting- 
time distribution seems to follow mainly a single ex- 
ponential distribution. It would therefore appear ap- 
pealing to parametrize this distribution and to draw 
the waiting times from the analytical expression of 
the parametrization. In this way the propagation of 
the auxiliary states could be avoided and a speedup of 
the propagation by a factor of two could be gained. 
However, it is not clear if the waiting times of the 
Kohn-Sham system follow always an exponential dis- 
tribution. In particular, the shape of the distribution 
is unknown, when, e.g., ionic motion is involved, or 
when the system is subject to strong external electric 
or magnetic fields. Therefore, to be on the safe side, 
in the present work we always sample the waiting- 
time distribution by propagating the auxiliary sys- 
tem. 



It is also worth noting that the average over 
stochastic realizations of the ensemble generally con- 
verges faster when the system-bath interaction in- 
creases. In the opposite limit the convergence is slow. 
When the system-bath interaction is very weak, only 
a small damping will be exerted by the S term in 



Eq. ( 56 1 , and hence it takes longer for the norm of the 



auxiliary wave functions to drop below their waiting 
times. This in turn implies that fewer jumps occur 
and hence more stochastic realizations are required 
to converge to a smooth observable distribution. 

4. Applications 

4-.1. Finite Systems 

In the last section we have introduced technical de- 
tails for the quantum jump algorithm that we use to 
simulate the stochastic process associated with the 
stochastic Kohn-Sham equation, Eq. (47). In this 



section we apply the algorithm to molecular systems 
with and without clamped ions. As first example 
we consider a situation of clamped ionic coordinates. 
As testcase we investigate a (l,4)-phenylene- linked 
zincbacteriochlorin-bacteriochlorin complex. Due to 
an extra Mg atom in the left porphyrin ring of the 
complex the molecule is not fully symmetric. As a re- 
sult, the highest-occupied molecular orbital (HOMO) 
of this molecule is located on the left porphyrin 
ring, whereas the lowest-unoccupied molecular or- 
bital (LUMO) is located on the right porphyrin ring, 
cf. Fig. [2j This system has been used as a model to 





Figure 2: Real part of the HOMO (left panel) 

and LUMO (right panel) orbitals of (l,4)-phenylene-linked 
zincbacteriochlorin-bacteriochlorin. 

study charge-transfer excitations in linear response 
TDDFT [32] ■ Here instead we consider open and 
closed system real-time propagation. We prepare the 
zincbacteriochlorin-bacteriochlorin complex in an en- 
tangled initial-state, where the orbital of the HOMO 
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Figure 3: Snapshots of the time evolution of the HOMO or- 
bital of zincbactcriochlorin-bacteriochlorin with clamped ions. 
The plots display the real part of the orbital. In the left col- 
umn the closed quantum system evolution is shown at different 
points in time and the right column displays the evolution of 
the system with a coupling to a thermal bath. A rather fast re- 
laxation rate of t = 1 a.u. has been used for the bath operator 
in the open quantum system case. 



is replaced by 

TDKS 
V>HOMO (t = 0) = 



1 

71 



■ GS 
^HOMO 



r / GS 



(56) 



where V'ho 



mo and "0LUMO denote the ground-state 
HOMO and LUMO, respectively. For all other or- 
bitals the ground-state configuration is used at the 
initial time. Starting from this excited initial Slater 
determinant the system is then evolved freely in time 
without any external fields. For the bath operators 



we employ the model of Eq. ( 48 1 introduced in section 
2.5| at zero temperature. The dynamics of the system 
is illustrated in Fig. [3] where we plot the real part 
of the HOMO orbital for different snapshots in time. 
The left panel summarizes the closed system evolu- 
tion and the right panel an ensemble average over 



100 stochastic realizations in the open system case. 
Let us first focus on the closed quantum system case. 
Due to the entangled initial state in Eq. (56), the 
time-dependence of the orbital V'homo (*) has mainly 
oscillatory phase contributions exp(— zchomo^) and 
exp(— ieLUMO*) from the ground-state HOMO and 
LUMO, respectively. Only small nonlinearities arise 
due to the dependence of the Kohn-Sham Hamilto- 
nian on the time-dependent density. Since the system 
is propagated as closed quantum system, the phase 
oscillations would continue indefinitely. 

On the other hand, the open auantum system evo- 
lution in the right panel of Fig. p shows for ip HO uo (t) 
a fast relaxation from the entangled initial state back 
to the HOMO which is localized on the left porphyrin 
ring. If we now imagine computing Ehrenfest forces 
from these orbital contributions, it is clear that the 
forces will differ qualitatively in the closed and open 
quantum system cases. While in the closed quan- 
tum system case the forces will be oscillatory, they 
will show relaxation behavior similar to the orbitals 
in the open quantum case. This example emphasizes 
that the coupling of electronic degrees of freedom to a 
thermal bath yields qualitatively different forces com- 
pared to standard QMD approaches. 

This observation motivates our second example, 
where we consider a stochastic QMD simulation for 
a neon dimer. In this case the ions are not clamped 
at the equilibrium configuration. Instead we use 
stretched initial positions for the ions of the dimer 
as initial state for the open and closed system propa- 
gation. If we would treat the ions quantum mechan- 
ically, then the bath operators would also act on the 
nuclear wavefunctions. However, since we have re- 
stricted ourselves here to the limit of classical ions 
we replace this action of the bath operators by mod- 
ifying the velocities of the ions. At every occasion 
when the bath operators act on the electronic wave- 
functions we draw new velocities for the ions from 
a Maxwcll-Boltzmann distribution. This is a simple 
approach but can be improved with e.g., recently in- 
troduced stochastic thermostats for the ions [22] . 

In the left panel of figure Fig. [4] we show the ionic 
positions of the dimer as function of time for a stan- 
dard Ehrenfest TDDFT closed quantum system evo- 
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Figure 4: Left panel: Here we show the ionic positions of a Neon dimer as function of time for a closed quantum system. As 
initial condition we have selected a stretched configuration of the dimer which results in an indefinite coherent oscillation of the 
two nuclei. Right panel: Using the same initial state we have evolved with SQMD a stochastic ensemble of trajectories. Shown 
is the average of the nuclear positions for an ensemble with 100 stochastic realizations. As relaxation time for the simulation 
we have employed r = 300 fs. The ionic velocities follow a Maxwell-Boltzmann distribution with a temperature of 290K. 



lution. In the right panel we display for the same ini- 
tial conditions the open quantum system evolution 
within SQMD. For the SQMD simulation we have 
employed a relaxation time of r = 300 fs and an av- 
erage over 100 stochastic trajectories has been per- 
formed which results in a smooth decay of the nuclear 
oscillations. 

4-2. Extended Systems 

So far we have considered only finite molecular sys- 
tems. However, a large class of applications requires 
also the treatment of periodic boundary conditions in 
one, two or three dimensions. This includes for ex- 
ample decoherence and dissipation in nanowires, elec- 
tronic relaxation on surfaces, or hot electron thermal- 
ization in bulk systems. For these cases it is desirable 
to extend our approach to periodic systems. In this 
section we briefly discuss the necessary steps in order 
to apply SQMD to extended systems. 

There are some extra details and conditions that 
have to be satisfied in order to treat periodic systems 
with SQMD. As a first step we expand the stochastic 
Kohn-Sham orbitals of the periodic system of interest 
in the complete set of the corresponding ground-state 
Bloch-orbitals ^fc(r) 



This gives rise to stochastic expansion coefficients 
dk(t) which are then propagated in time using, e.g., 
the quantum jump algorithm that we have presented 
in section 3.2 Similar to the case of molecules, the 



ipk(L,t) = ^<2 fe '( ^ )¥ : >fe'(r)• 



(57) 



fermionic nature of the electronic subsystem needs to 
be taken into account by orthogonalizing the occu- 
pied states after each application of the bath operator 
(cf. step (4) of the quantum jump algorithm). 

In addition, care has to be taken for the choice 
of gauge for the vector and scalar potentials in the 
Hamiltonian of the extended system. Here, the 
same restrictions apply as in standard closed-system 
TDDFT simulations. In practice, we consider only 
purely time-dependent vector potentials which retain 
the periodicity of the considered system at all times. 
In the present context we have to assume in addition 
that the bath operators retain the periodicity of the 
extended system as well. This restricts the choice of 
baths represented by local operators that satisfy the 
condition 

S(r) = S(r + R) (58) 

where R denotes the usual displacement vector of the 
unit cell. This may exclude certain relaxation mech- 
anisms. However, the importance of these relaxation 
and dephasing channels can always be checked by in- 
creasing the size of the supercell that is used in the 
simulation. 
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While we do not have fully implemented this 
scheme yet, we want to argue about important physi- 
cal processes that can be studied with this approach. 
For instance, one could study adsorption of molecules 
on surfaces whose opposite side is set on a thermal 
stage that keeps the electron and/or ion temperatures 
fixed. Again, this could be accomplished in a super- 
cell geometry by coupling some "bulk" layers away 
from the surface with a local operator that maintains 
energy equilibrium in that region (an example of such 
operator is given in Ref. [13]). The rest of the sys- 
tem is let to follow its own dynamics. If we excite 
the molecules and/or surface - e.g., by application of 
a short electromagnetic field - electrons and ions can 
then distribute energy and momentum first in the lay- 
ers adjacent to the surface and then relax energy into 
the bath, where they would thermalize to the appro- 
priate canonical distributions. Analogously, we could 
monitor energy relaxation of electrons and ions in a 
bulk exited either thermally or electrically and kept 
at a given temperature by a thermal stage. Impor- 
tant phenomena that are then accessible would be, 
e.g., phase transitions driven by dissipative effects. 

5. Conclusions 

In summary, we have presented a detailed account 
of stochastic quantum molecular dynamics. The ap- 
proach is based on a stochastic Schrodinger equation, 
which may or may not describe Markovian dynam- 
ics - although we have focused the discussion to the 
Markovian case. Our approach allows us to describe 
the dynamics of electrons and ions coupled to one 
or many external environments. For simplicity we 
have restricted our examples to the situation of clas- 
sical ions, but the approach is, in principle, valid also 
for quantum ions. Although we have not reported 
any actual implementation of SQMD for periodic sys- 
tems, we have outlined the theory behind its exten- 
sion to this important case. Work along these lines 
is in progress and will be reported elsewhere 14J. 

This approach is thus amenable to studying many 
interesting phenomena related to energy relaxation 
and dephasing of the electronic subsystem in the pres- 
ence of ionic dynamics, such as local ionic and elec- 
tronic heating in laser fields, relaxation processes in 



photochemistry, etc., a feature that is lacking in any 
"standard" molecular dynamics approach. 

We acknowledge support from DOE under grant 
DE-FG02-05ER46204 and Lockheed Martin. 
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